function syns = nn_train(X, Y, nn_size)
%{
Sample 
Input:
    X = [0 0 1; 0 1 1; 1 0 1; 1 1 1];
    Y = [0 1 1 0]';
    nn_size = [3 4 1];

syns{1} =
    5.0983    3.2440   -2.2129    7.6981
    5.2360   -6.3985    4.7233   -5.6545
   -0.8552   -0.6890    0.4585    2.1030
syns{2} =
    9.0207
    9.1934
   -3.3483
  -11.3997
%}
    train_times = 60000;

    %% varify input arguments
    layers_len = length(nn_size);
    assert(nn_size(1) == size(X,2));
    assert(nn_size(end) == size(Y,2));

    %% init syns
    syns = cell(1, layers_len - 1);
    for i = 1:layers_len-1
        syns{i} = 2 * rand(nn_size(i), nn_size(i+1)) - 1;
    end
    
    %% train
    layers = cell(1, layers_len);
    layers_error = cell(1, layers_len);
    layers_delta = cell(1, layers_len);
    for i = 1:train_times
        % forward
        layers{1} = X;
        for j = 2:layers_len
            layers{j} = sigmoid(layers{j-1} * syns{j-1});
        end
        
        % backward
        layers_error{end} = Y - layers{end};
        layers_delta{end} = layers_error{end} .* sigmoid(layers{end}, true);
        for j = layers_len-1:-1:1
            layers_error{j} = layers_delta{j+1} * syns{j}';
            layers_delta{j} = layers_error{j} .* sigmoid(layers{j}, true);
        end
        
        % update syns
        for j = layers_len-1:-1:1
            syns{j} = syns{j} + layers{j}' * layers_delta{j+1};
        end
    end
end
